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ABSTRACT 

Self diffusion coefficients and binary Maxwell-Stefan diffusion coefficients were determined 
by equilibrium molecular dynamics simulations with the Green-Kubo method. The study 
covers five pure fluids: neon, argon, krypton, xenon, and methane and three binary mix- 
tures: argon-|-krypton, argon-|-xenon, and krypton-|-xenon. The fluids are modeled by 
spherical Lennard-Jones pair-potentials, with parameters which were determined solely 
on the basis of vapor-liquid equilibria data. The predictions of the self diffusion coeffi- 
cients agree within 5% for gas state points and about 10% for hquid state points. The 
Maxwell-Stefan diffusion coefficients are predicted within 10%. A test of Darken's model 
shows good agreement. 
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1. INTRODUCTION 

Diffusion plays an important role in many chemical processes, such as catalysis or adsorp- 
tion. On the other hand, the measurement of diffusion coefficients is a time consuming 
and difficult task [1]. Molecular simulation offers the possibility to straightforwardly de- 
termine diffusion coefficients on the basis of a given molecular model. Both self diffusion 
coefficients and binary Maxwell-Stefan (MS) diffusion coefficients can be obtained by non- 
equilibrium molecular dynamics (NEMD) or equilibrium molecular dynamics (EMD). In 
this work, EMD is chosen. 

From the pioneering work of Alder and Wainwright with hard spheres [2]|3] . the simula- 
tion of diffusion coefficients has been an area of continuous research. There are several 
contributions in which self diffusion coefficients binary [71IS|U][Tmillll2j and ternary 

diffusion coefficients [T3p^ for noble gases, methane and mixtures of these are calculated. 
Less frequently, investigations with multi-center Lennard- Jones models, e.g. mixtures of 
CH4+SF6 [15J and CH4+CF4 [16j, or polar ffuids [T7IITH] have been performed. With 
the exception of Refs. 5 and 6, all investigations from the literature cover only diffusion 
coefficients in the liquid phase and only for a limited range of state points. 

This is the aim of the present work in which, as a first step, only simple ffuids are con- 
sidered. Self diffusion coefficients for five pure ffuids: neon, argon, krypton, xenon, and 
methane (both liquid and gas) and three binary mixtures: argon-|-krypton, argon+xenon, 
and krypton+xenon (gas) are predicted based on molecular models from the literature and 
compared with experimental data. The pure component parameters of these models were 
determined from vapor-liquid equilibria data alone [19] . Binary mixtures were modeled us- 
ing one adjustable parameter for the unlike interaction which was fitted to vapor-pressure 
data of the mixtures [20121] • Throughout the present work, for the molar mass the stan- 
dard value from the literature is used [22] ■ The simulation results on diffusion coefficients 
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from the present work are therefore predicted from vapor-hquid equihbria alone and ob- 
tained without any adjustment to diffusion or other transport data. The studied systems 
are those for which both molecular models and experimental data were available. 
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2. METHOD 

2.1. Molecular Models 

In this work, only noble gases and methane are considered, so that the description of the 
molecular interactions by the Lennard- Jones 12-6 (L J) potential is sufficient and physically 
meaningful. The LJ potential u is defined by 

Uij{r) = 4e,j 

where i and j are the species indices, aij is the LJ size parameter, eij the LJ energy 
parameter, and r the center-center distance between two molecules. Pure substance pa- 
rameters (jjj and en are taken from Ref. 19 as given in Table [B They were adjusted by 
Vrabec et al. [I9l to experimental pure substance vapor-liquid equilibrium data alone. For 
modeling mixtures, parameters for the unlike interactions are needed. Following previous 
work [20ll2T] . they are given by a modified Lorentz-Berthelot combination rule 




(q-ii + 0-22) 

f^i2 = ^ , (2) 

and 

ei2 = ^ ■ y/eue^. (3) 
where ^ is an adjustable binary interaction parameter. This parameter allows an accurate 
description of the binary mixture data [20|2T] and was determined by an adjustment to 
one experimental vapor-liquid equilibrium state point. As for the mixture krypton+xenon 
no binary interaction parameter is available, ^ =1 was assumed. The binary interaction 
parameters are listed in Table [Tll 
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2.2. Diffusion Coefficients 

Diffusion coefficients can be calculated by equilibrium molecular dynamics through the 
Green-Kubo formalism [2^11^ . In this formalism, transport coefficients are related to 
integrals of time- correlation functions. There are various methods to relate transport 
coefficients to time-correlation functions; a good review can be found in [25]. The self 
diffusion coefficient Di is given by [13] 



where vf(t) expresses the velocity vector of molecule k of species i and the notation 
< ... > denotes the ensemble average. Equation (jH) yields the self diffusion coefficient for 
component i averaging over Ni molecules. The expression for the binary Maxwell- Stefan 
diffusion coefficient D12 is given by [13] 



where Mi denotes the molar mass of molecules of species A'^i the number of molecules 
of species 1 and xi and X2 are the mole fractions. 

To compare MS diffusion coefficients to available experimental data, it is necessary to 
transform the MS diffusion coefficients to Pick diffusion coefficients. There is a direct 
relation between binary MS diffusion coefficients D12, and binary Pick diffusion coefficients 
D12 [2S], which is given by 




(4) 




D12 = Du ■ g, 



(6) 



with 




(7) 



7 



where /zi is the chemical potential of species 1, ks the Boltzmann constant and T the 
temperature. 

Because the present simulations provide both binary MS diffusion coefficients and self 
diffusion coefficients, it is possible to test the often used model of Darken [27f28] . It gives 
an estimate of the MS diffusion coefficient, D^2i from the self diffusion coefficients of both 
components in a binary mixture Di and D2, 

D% = Di-xi + D2- X2. (8) 

2.3. Simulation Details 

The molecular simulations were performed in a cubic box of volume V containing stan- 
dard = 500 molecules modeled by the LJ potential. The cut-off radius was set to 
Tc = 5cr; standard techniques for periodic boundary conditions and long-range corrections 
were used [22] • The simulations were started with the molecules in a face-center-cubic 
lattice with random velocities, the total momentum of the system was set to zero, and 
Newton's equations of motion were solved with a velocity- Verlet algorithm [29]. The time 
step for this algorithm was set to At • ^Je^Jraxj o\ = 0.001 for liquid and 0.01 for gas state 
points. The diffusion coefficients were calculated in a NVE ensemble, using Eqs. (jlj) and 
(E]). The relative fluctuation in the total energy in the NVE ensemble was less than 10^^ 
for the longest run, which yields a temperature drift of less than 0.5 K. The simulations 
are initiated in a NVT ensemble until equilibrium at the desired density and temperature 
is reached. 25 000 time steps were used for that equilibration. Once the equilibrium is 
reached, the thermostat is turned off, and then the NVE ensemble is invoked. The exper- 
imental data which was used for comparison to our simulations is often reported at given 
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pressure and temperature. In these cases, a prior isobaric-isothermal NpT simulation [30] 
was performed, from which the corresponding densities for the NVE ensemble were taken. 
The statistical uncertainty of the diffusion coefficients was estimated with the standard 
block average technique [31] . 

The self diffusion coefficient is a property related to one molecule. It is possible to obtain 
very good statistics with a few independent velocity autocorrelation functions (VACF). 
The self diffusion coefficients were calculated by averaging over 200 independent VACF 
each with 500 molecules, i.e. a total of 10^ VACF. For gas densities, the VACF decays very 
slowly and therefore long simulation runs were necessary in order to achieve the VACF 
decay and hence independent time-origins. Here, a compromise between simulation time 
and time-origin independence had to be made. In order to keep the simulation time low, 
and following the work of Schoen and Hoheisel [^, the separation between time origins 
was chosen at least as long as the VACF needs to decay to 1/e of its normalized value. The 
choice of this separation time and the production phase depended upon the temperature 
and density conditions. In theory, as Eq. (jlj) shows, the value of the diffusion coefficient is 
determined by an infinite time integral. In fact, however, the integral is evaluated based 
on the length of the simulation. The integration must be stopped at some finite time, 
ensuring that the contribution of the long-time tail [3J is small compared to the desired 
statistical uncertainty of the diffusion coefficient. Figure [1] shows the behavior of the 
VACF and its integral given by Eq. (jl]) for two selected gas state points of argon. As can 
be seen, for the higher density state point, the VACF has decayed after 500 ps to less 
than 1% of its normalized value. Later it oscillates around zero. The same can be seen 
after 1500 ps at the lower density state point. It was assumed here that the VACF has 
fully decayed when these oscillations reached a threshold below 0.5% of their normalized 
value. Furthermore, a graphical inspection of the VACF integral was made, in order to 
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verify a sufficient integration time. 

An important time scale to calculate the VACF is the time that a sound wave takes to cross 
the simulation box. VACF calculated for times higher than that may show a systematic 
error [321133] . That criterion was verified using the experimental speed of sound taken from 
[3^ . For the simulations of gases, the VACF decay time was found to be higher than that 
time. To check whether this to leads to a systematic error in the present simulations, 
the system size was varied. For the lowest density state point of argon, where the above 
mentioned problem would be expected to be most severe, simulations with a constant 
number of time origins and increasing system sizes were carried out. System sizes of 
= 864, 2048, 4000, and 6912 molecules were investigated. All results were found to 
agree within the statistical uncertainty, and no size dependence could be observed. It is 
therefore concluded that no systematic error due to system size in gas phase simulations 
is present. 

The Maxwell- Stefan diffusion coefficient is a collective quantity and therefore the statis- 
tics can only be improved by averaging over longer simulation runs. The MS diffusion 
coefficients were calculated by averaging over 2000 velocity correlation functions (VCF) 
as proposed by Schoen and Hoheisel [8]. In order to obtain independent time origins, 
similar criteria as employed for the self diffusion coefficients were used to determine the 
necessary length of the VCF. 
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3. RESULTS 

3.1. Self Diffusion Coefficients 

Figure [2] shows the results for the self diffusion coefficients of neon, argon, and krypton 
compared with experimental data for gas state points [3l]. The lines in Fig. [2] are the 
results of the correlation of Liu et al. [33] using the LJ parameters from Table [B The data 
are given at constant pressure and at different temperatures. Figure O shows the results 
for the self diffusion coefficients for neon, argon, krypton, xenon, and methane compared 
with experimental data [36|37f38f39f40j for liquid state points. The lines in Fig. [3] are the 
results of the correlation of Liu et al. [35j using the LJ parameters from Table [11 The 
data are given at constant temperature and at different pressures. Overall, a very good 
agreement between the predictions and the experimental data is found. The best results 
are found for neon, argon, and krypton in the gas phase with deviations within a few 
percent. The results for liquid state points show somewhat higher relative deviations from 
the experimental data (around 10%). It can be seen that the correlation agrees reasonably 
well with the simulation data, typical deviations are about 5%. This accuracy lies in the 
range claimed by the authors of Ref. 35. 

3.2. Binary Maxwell-Stefan Diffusion Coefficients 

Binary MS diffusion coefficients were calculated for the gaseous mixtures argon+krypton, 
argon+xenon, and krypton+xenon. The results are compared to experimental Fick diffu- 
sion coefficients. The thermodynamic factor Q, that relates the MS diffusion coefficient 
to the Fick diffusion coefficient, cf. Eq. ([6]), is assumed to be unity for all cases studied 
here. This is supported by the calculations of several authors [7]l8|llf41j . As a test, Q 
was estimated by two simulations to calculate a finite difference [12] for each mixture at 
the most dense state point. The assumption Q =1 was confirmed within the statistical 
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uncertainty of the calculations. 

Figure H] shows the simulation results for the mixture argon+krypton in comparison to 
experimental data taken from Ref. 43. The continuous line in Fig. H] are the results of 
the correlation of Darken [27f28] . In this case, the experimental data [44j were reported 
at constant temperature. Figure [5] shows the results for the mixtures argon+xenon and 
krypton+xenon at constant pressure. Good agreement between the predictions and the ex- 
perimental values is found. The best results are observed for the mixture argon+krypton. 
This mixture shows typical relative deviations of 4% from the experimental data; the 
corresponding numbers are 8% for argon+xenon, and 16% for krypton+xenon. It must 
be pointed out that for the mixture krypton+xenon, no binary interaction parameter ^ 
was available. 

In Figs, m and [5], it is interesting to analyze the performance of the empirical model of 
Darken for estimating MS diffusion coefficients in the gas phase on the basis of known 
binary self diffusion coefficients. The figures show that the model of Darken agrees very 
well with the binary MS diffusion coefficients, typical deviations are about 5%. 
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4. CONCLUSION 

In the present work, molecular models of simple fluids that were adjusted to vapor-liquid 
equilibrium only were used to predict self and MS diffusion coefficients. The diffusion 
coefficients were determined with molecular dynamics simulations using the Green-Kubo 
method. Five pure fluids and three binary mixtures were studied covering a broad range 
of state points. The fluids were modeled with the Lennard- Jones pair potential with 
parameters taken from the literature. It is found that the prediction of diffusion coefficients 
from vapor-liquid equilibrium data using that simple model yields good results. This 
supports the finding that the spherical LJ 12-6 potential is an adequate description for 
the regarded noble gases and methane. When molecular models are adjusted to diffusion 
coefficient data, excellent descriptions can be expected. It is worthwhile to extend the 
study to more complex fluids. 
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Table I 

Potential Model Parameters for the Pure Fluids used in this Work"'^. 



Fluid 






M (p- ■ mn\~^\ 


Neon 


o cm n 




on 1 on 


Argon 


o.oyoz 


lib. i\) 


on n/i c 


Krypton 




IDZ.Oo 


oo.o 


Xenon 


3.9011 


227.55 


131.29 


Methane 


3.7281 


148.55 


16.043 



^ Values taken from Ref. 19. 

^ The Molar Mass M was taken from Ref. 22. 
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Table II 

Binary Interaction Parameters for the Binary Mixtures taken from Ref. 21. 



Mixture 




Argon + Krypton 


0.988 


Argon + Xenon 


1.000 


Krypton + Xenon 


0.989 
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List of Figures 



1 Large plot: Velocity autocorrelation functions. Small plot: Integral 
following Eq. (jH). Both plots are shown for selected state points of argon: 

T=77.7 K and p=163.26 mol ■ m'^; T=353.2 K and p=34.49 

mol ■ m~^. 

2 Self diffusion coefficients of neon, argon and krypton (gas phase) predicted 
by molecular simulation compared to experimental data [M] at j9=0.1013 
MPa. neon: □ exp., ■ sim.; argon: A exp., A sim.; krypton: o exp., • sim. 
The solid lines represent the results of the correlation of Liu et al. 

3 Self diffusion coefficients of neon, argon, krypton, xenon, and methane 
(liquid and gas phase) predicted by molecular simulation compared to 
experimental data [36|37f38f39f40j at constant temperatures and different 
pressures, neon T=37 K: □ exp., ■ sim.; argon T=323 K: A exp., ▲ 
sim.; krypton T=273 K: o exp., • sim.; xenon T=298 K: () exp., ^ sim.; 
methane T=298 K: V exp., ▼ sim. The solid lines represent the results of 
the correlation of Liu et at |35] . 

4 Binary Maxwell-Stefan diffusion coefficients for gaseous equimolar 
mixtures of argon+krypton predicted by molecular simulation compared 
to experimental data [13] at T=323.16 K : argon+krypton A exp., ▲ sim. 
The solid lines represent the results of the correlation of Darken [.27|28j. 
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Binary Maxwell-Stefan diffusion coefficients for gaseous equimolar mixtures 
of argon+xenon and krypton+xenon predicted by molecular simulation 
compared to experimental data at ]9=0.1013 MPa: argon+xenon: o exp., 
• sim.; krypton+xenon V exp., T sim. The solid lines represent the results 
of the correlation of Darken [27I2H]- 
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Fig. 1. 
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Fig. 2. 
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Fig. 3. 




Fig. 4. 
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Fig. 5. 
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